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Invited  Paper 


Abstract — Realistic  and  accurate  numerical  simulations  of  elec¬ 
trostimulation  of  tissues  and  full-body  biomodels  have  been  de¬ 
veloped  and  implemented.  Typically,  whole-body  systems  are  very 
complex  and  consist  of  a  multitude  of  tissues,  organs,  and  subcom¬ 
ponents  with  diverse  properties.  From  an  electrical  standpoint, 
these  can  be  characterized  in  terms  of  separate  conductivities  and 
permittivities.  Accuracy  demands  good  spatial  resolution;  thus, 
the  overall  tissue/animal  models  need  to  be  discretized  into  a  fine¬ 
grained  mesh.  This  can  lead  to  a  large  number  of  grid  points 
(especially  for  a  three-dimensional  entity)  and  can  place  pro¬ 
hibitive  requirements  of  memorv  storage  and  execution  times  on 
computing  machines.  Here,  the  authors  include  a  simple  yet  fast 
and  efficient  numerical  implementation.  It  is  based  on  LU  de¬ 
composition  for  execution  on  a  cluster  of  computers  running  in 
parallel  with  distributed  storage  of  the  data  in  a  sparse  format. 
In  this  paper,  the  details  of  electrical  tissue  representation,  the 
fast  algorithm,  the  relevant  biomodels,  and  spedfic  applications 
to  whole-animal  studies  of  electrostimulation  are  discussed. 

Index  Terms — Distributed  currents,  LU  decomposition,  parallel 
computing,  tissue  modeling,  whole  body. 

I.  Introduction 

Electrical  excitation,  which  has  been  used  to  stimulate 
both  the  central  and  the  peripheral  nervous  system,  has 
a  variety  of  potential  diagnostic  and  therapeutic  applications 
[l]-[3].  For  example,  electrically  stimulated  neurogenesis  is  a 
potential  tool  for  enabling  the  production  of  new  nerve  cells 
from  neuronal  stem  cells  [4],  [5].  It  is  used  in  implantable  de¬ 
vices  for  neuromuscular  stimulation;  these  devices  are  designed 
to  control  the  contraction  of  paralyzed  skeletal  muscles,  thereby 

Manuscript  received  January  7,  2006;  revised  March  31,  2006.  This  work 
was  supported  in  part  by  the  National  Training  and  Information  Center  (NTIC) 
and  in  part  by  the  Air  Force  Office  of  Scientific  Research  (AFOSR),  U.S.  Air 
Force. 

A.  Mishra,  R.  P.  Joshi,  and  K.  H.  Schoenbach  are  with  the  Department  of 
Electrical  and  Computer  Engineering,  Old  Dominion  University.  Norfolk.  VA 
23529-0246  USA  (e-mail:  rjoshi(§>odu.edu). 

C.  D.  Clark.  HI.  is  with  General  Dynamics.  Radio  Frequency  Radiation 
Branch.  AFRL/HEDR.  Brooks  City  Base,  San  Antonio,  TX  78235-5417  USA 
(e-mail:  Clifton.Clark.ctr@ brooks.af.mil). 

Color  versions  of  Figs.  I  and  3-19  are  available  online  at  http:// 
ieeexplore.org. 
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producing  functional  movements  in  patients  with  stroke  and 
spinal  cord  injuries  [6],  [7].  Electrical  excitation  is  al.so  a  u.seful 
tool  to  study  the  properties  and  functions  of  nerves  (including 
the  brain)  and  muscles  (including  the  heart).  It  also  provides 
information  on  strength-duration  characteristics,  including  the 
inherent  time  constants.  Excitation  with  ultrashort  pulses  is 
also  an  important  issue  for  the  health  and  safety  assessment 
of  ultrawideband  (UWB)  sources,  which  produce  nanosecond 
pulses.  This  issue  has  been  discussed  at  length  elsewhere  [8]. 

In  general,  muscle  excitation  can  be  achieved  either  remotely 
through  the  principle  of  electromagnetic  induction  [9]-[  1 1  ] 
or  directly  through  electrical  contact  [12].  UWB  pulses  have 
not,  to  date,  shown  significant,  robust,  and  reliable  biological 
effects.  Such  UWB  studies  using  rats  probed  behavioral  tera¬ 
tology,  heart  rate,  blood  pressure,  brain  histology,  and  genetic 
alterations  [13]-[16].  This  contribution  focuses  on  pulsed  exci¬ 
tation  delivered  to  biological  tissues  (and  whole  bodies)  through 
direct  electrical  contacts.  Thus,  it  is  assumed  that  electrodes 
can  be  directly  applied  on  the  muscle/tissue  surface,  and  motor 
nerve  fibers  within  the  muscle  are  then  excited  by  the  potential 
created  within  the  muscle  by  the  external  source.  In  general, 
the  potential  can  also  be  applied  through  a  conductive  medium 
surrounding  the  biomass  as  discussed  in  a  previous  report  [17]. 

The  use  of  ultrashort  electrical  pulses  in  this  context  is  an 
emerging  topic  of  interest  [  1 8].  Such  pulses  of  nanosecond  du¬ 
ration  have  been  shown  to  have  the  ability  to  penetrate  the  outer 
(plasma)  membrane — without  affecting  it — to  create  large 
trans-membrane  potentials  across  subcellular  organelles  |19]. 
Thus,  for  example,  triggering  of  neurotransmitter  or  calcium 
release  is  possible  through  the  use  of  ultrashort  pulses  [20].  Use 
of  this  new  technology  to  study  submicrosecond  pulse  widths 
might  reveal  new  biological  phenomena. 

The  current  study  was  driven,  in  part,  from  health  and  safety 
concerns  of  electrostimulation.  The  overall  goal  is  to  under¬ 
stand,  quantify,  and  predict  possible  adverse  electrophysiolog- 
ical  changes,  including  electrically  induced  organ  failures  in 
humans  and  whole  animals,  due  to  electrostimulation.  As  first 
di.scussed  by  McNeal  [21  ],  this  overall  objective  requires  a  two- 
step  approach.  The  first  essential  component  is  the  development 
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Fig.  1.  Rai  rmxiel  enclosed  in  a  simulation  box.  Dimensions  arc  in 
millimeters. 


of  a  numerical  model  for  quantifying  the  microscopic  currents 
and  electrically  induced  potentials  due  to  an  external  voltage 
pulse  in  whole-body  systems.  Such  capability  can  also  provide 
a  useful  tool  for  optimal  electrode  design  and  placement.  The 
second  step  is  to  use  the  derived  excitation  potentials  within 
the  whole-body  system  to  determine  the  biological  response  of 
nerves,  muscles,  and  related  electrochemical  reactions.  Here,  in 
this  contribution,  our  focus  is  on  the  first  part — the  development 
of  an  efficient  numerical  model  for  predicting  the  voltage- 
induced  currents  and  potentials. 

The  whole-body  system  is  very  complex  and  consists  of  a 
multitude  of  tissues,  organs,  and  subcomponents  with  diverse 
properties.  From  an  electrical  standpoint,  these  can  all  be 
characterized  in  terms  of  separate  electrical  properties,  such  as 
conductivity  and  permittivity.  Accuracy  demands  good  spatial 
resolution,  and  so  the  overall  whole-animal  system  needs  to  be 
discretized  into  a  fine-grained  mesh  of  individual  grid  points. 
The  result,  especially  for  a  three-dimensional  (3-D)  entity,  is 
that  a  very  large  number  of  grid  points  would  typically  be 
required  for  the  numerical  computations.  The  requirements  for 
memory  storage  and  execution  times  on  computing  machines 
can  quickly  become  prohibitive  in  such  situations.  Here,  we 
include  a  simple  yet  fast  numerical  implementation  for  ex¬ 
ecution  in  parallel  on  a  system  of  computers  (parallelized 
cluster  computing).  The  details  of  this  fast  algorithm  and  its 
implementation  are  given,  followed  by  a  specific  application  to 
whole-animal  studies  of  electrostimulation. 

II.  Simulation  Method 

In  the  present  approach  for  voltage  and  current  calculations, 
the  entire  simulation  volume  is  first  discretized  into  a  dis¬ 
tributed  array  of  uniformly  spaced  grid  points  within  a  cubic 
simulation  volume.  The  simulation  volume  encompasses  the 
entire  bioma.ss  (i.e.,  full-animal  body)  under  study.  Due  to  the 
nonrectangular  and  curved  nature  of  the  full  animal  objects 
placed  within  the  simulation  box,  regions  of  surrounding  air 
are  also  included  in  the  cubic  region.  Thus,  irrespective  of 
the  actual  irregular  shape  of  a  body,  each  raw  numerical  file 
involves  a  regular  encompassing  volume,  as  illustrated  in  Fig.  1 . 
This  figure  shows  a  volumetric  visualization  of  the  entire  rat 


Fig.  2.  Schematic  of  the  di.scrciized  electrical  model  used. 

data  present  in  the  corresponding  raw  model  file.  Fig.  1  also 
highlights  different  tissue  types  by  different  colors.  The  en¬ 
tire  simulation  box  is  divided  into  a  Cartesian  grid  of  points 
(Nx*  Ny,  and  Nz  along  the  three  axes).  Thus,  each  grid  point 
has  at  most  a  total  of  six  nearest  neighbors.  For  electrical  com¬ 
putations,  the  nearest  neighbor  sub-regions  were  represented 
in  terms  of  a  parallel  resistor-capacitor  combination  to  account 
for  both  conduction  and  displacement  current  flows  within  the 
biosystem.  In  essence,  a  distributed  RC  model  interconnected 
to  (at  most)  six  neighbors  was  used  for  the  electrical  repre¬ 
sentation.  This  is  shown  schematically  in  Fig.  2.  The  applied 
voltage  was  taken  to  be  a  user-specified  time-dependent  input 
waveform  for  the  numerical  solution,  and  we  essentially  im¬ 
posed  Dirichlet  boundary  conditions  locally.  Neumann  bound¬ 
ary  conditions  of  zero  net  current  flow  were  applied  on  all 
other  boundary  regions  that  did  not  contain  any  electrodes.  For 
simplicity,  the  passive  RC  elements  were  chosen  to  be  linear, 
and  hence,  had  fixed  values  based  on  the  conductivity  [a)  and 
permittivity  (sr)  characteristics  of  the  local  region.  Inasmuch  as 
the  entire  body  consists  of  a  large  multitude  of  tissue  types, 
a  range  of  cr  and  e  values  was  obtained  from  actual  data* 
and  stored  in  a  look-up  table.  As  an  example.  Table  1  gives 
the  electrical  properties  of  some  tissues  that  were  used  in  this 
modeling  study. 

The  application  of  current  continuity  at  each  node  then 
results  in  a  set  of  coupled  equations  for  the  node  voltages, 
which  can  be  solved  at  each  time  step.  Time-dependent  values 
of  the  potential  across  each  discretized  subregion  and  current 
distributions  can  then  be  directly  obtained.  In  theory,  the 
voltages  can  be  obtained  through  Kirchhoff's  node  analysis.  At 
each  node,  equations  of  the  type 

E  [{Ae/L)d{^V]/dt  +  { AK}(Ac7/L)]  =0  ( 1 ) 

apply  with  A  being  the  elemental  cross-sectional  area  and  L  the 
length  of  each  segment.  These  can  be  cast  into  the  matrix  form 

[M]  [AKie+rf,  -  AV^I,]  =  [B{t)]  (2) 

‘  [Online].  Available:  flp://stan’iew.brooks.af.mil/EMF/do.simeir\ .model.s/ 
computer_coded_binary_files/ 
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TABLE  I 

Electrical  Parameters  for  Various  Constituent  Tissues 


TissucNamc 

Tissue  ID 

Sigma 

Epsilon  (real 
part) 

AlR.(c.\tcmaJ) 

0 

0 

I.OOE+OO 

AIR.  (internal) 

1 

0 

l.OOE+OO 

BILE 

2 

1.4 

4.00E^00 

BODY.FLUID 

3 

1.5 

4.00E+00 

EYE.  (cornea) 

4 

0.4 

4.00E+00 

FAT 

5 

0.01 

2.50E+00 

LYMPH 

6 

0.5 

4.00E+00 

MUSCOUS.MEMBRANE 

7 

0.0004 

4.00E-b00 

N  A I L  .S .  ( toes .  and .  fi nge rs ) 

8 

0.02 

2.50E+00 

NERVE. (spine) 

11 

0.006 

4.00E-K)0 

MUSCLE. 

17 

0.2 

4.00E+00 

HEART 

25 

0.05 

4.00E+00 

WHITE.MATTER 

30 

0.02 

4.00E+00 

STOMACH 

48 

0.5 

4.00E4-00 

GLANDS 

49 

0.5 

4.00E+00 

BLOODVESSEL 

65 

0.25 

4.00E+00 

LIVER 

68 

0.02 

4.00E+00 

GALL  BLADDER 

88 

0.9 

4.00E+00 

SPLEEN 

108 

0.03 

4.00E+00 

CEREBELLUM 

no 

0.04 

4  00F+00 

BONE. (cortical) 

111 

0.02 

2 50L+00 

CARTILAGE 

133 

0.15 

4.00E+00 

LIGAMENTS 

142 

0.25 

4.00E-f00 

SKIN/DERMLS 

143 

0.0002 

4.00E-^00 

INTESTINE.(largc) 

148 

0.01 

4.00E4-00 

T(X)TH 

152 

0.02 

2.50L+00 

GRAY.MATTER 

160 

0.02 

4.00E+00 

EYE.(lens) 

163 

0.3 

4.00E-H)0 

LUNG. (outer) 

164 

0.2 

4.00E+00 

INTESTINE.(small) 

168 

0.5 

4.00E+00 

EYE.(sclera/wall) 

183 

0.5 

4.00E4-00 

LUNG. (inner) 

184 

0.03 

2.50E+00 

PANCREAS 

188 

0.5 

4.00E{-H)0 

BLOOD 

189 

0.7 

4.00E+00 

CEREBRALSPINALFLUID 

190 

2 

4.00t+00 

EYE. (retina) 

203 

0.5 

4.00E+00 

EYE. (aqueous. humour) 

204 

1.5 

4.00E+00 

KIDNEYS 

207 

005 

4  00E+00 

BONE.MARROW 

209 

0.0005 

2.50E+00 

BLADDER 

227 

0.2 

2.50E+00 

TESTICLES 

228 

0.4 

4.00E+00 

PERFECT.CONDUCTOR 

229 

-1 

-I.OOE+OO 

2/3.MUSCLE. 

249 

0.1333 

2.70E+00 

PVC 

250 

0 

2.46 

FOAM 

251 

0.0004 

1.16 

TEM.(oId) 

252 

1.16 

1.95 

BONE,  (cancellous) 

253 

0.07 

2.50E+00 

TEM.(new) 

254 

1.3 

5.6 

UTAH 

255 

0.859 

39.2 

where  [M]  is  a  matrix  (asymmetric  in  general),  whereas  AV\t 
and  are  column  vectors.  Inasmuch  as  each  node  is 

connected  only  to  six  (at  most)  neighbors,  the  matrix  [M] 
is  extremely  sparse.  The  asymmetry  in  electrode  placement 
makes  [M]  nonsymmetric  in  location  and  values.  Fig.  3 
highlights  this  sparsity  by  way  of  a  plot  of  the  total  elements 
with  respect  to  the  nonzero  entries  for  cubic  volumes  of  varying 
size.  For  efficient  data  storage,  one  does  not  really  need  to  set 
up  the  matrix  in  its  entirety,  but  instead  only  record  the  small 
number  of  nonzero  elements. 

Mathematically,  the  desired  solution  at  each  time  step 
can  be  obtained  by  inverting  [M],  and  this  needs 


Fig.  3.  Sparsity  of  the  resulting  coefficient  matrix  [A/].  Nonzero  values 
(— X— )  (in  blue)  and  total  elements  (— *— )  (including  zeix)s)  for  different 
cases. 


to  be  carried  out  only  once  at  the  beginning.  However,  from 
a  practical  standpoint,  if  the  total  number  of  nodes  becomes 
large  (for  fine  resolution  within  a  whole-body  system),  then 
the  inversion  process  becomes  prohibitively  time  consuming 
and  eventually  intractable.  A  practical  solution  is  to  make  use 
of  the  LU  decomposition  of  a  matrix.  For  further  speed  up  and 
to  provide  for  the  large  memory  storage  necessary,  the  optimal 
solution  is  to  implement  the  LU  decomposition  on  a  cluster  of 
computers  running  in  parallel  with  distributed  storage  of  the 
data  in  a  sparse  format.  Computational  schemes  of  such  nature 
have  recently  been  devised  [221-[29],  and  here  we  discuss  our 
numerical  implementation. 

The  nonsymmetric  nature  of  the  matrix  [M]  eliminates 
the  possibility  of  using  the  Cholesky  factorization  method 
[30],  which  is  known  to  be  one  of  the  most  efficient  and  fast 
techniques.  The  extremely  sparse  structure  of  the  coefficient 
matrix  [M]  led  us  to  adopt  SuperLU.DIST.  This  is  a  scalable 
distributed-memory  sparse  direct  solver  for  asymmetric  linear 
systems  [29]  for  obtaining  the  node  voltages.  The  details 
of  the  algorithm  can  be  found  elsewhere  [29].  This  is  a 
numerically  stable  approach,  and  both  scalable  implementation 
and  scheduling  optimization  are  possible  because  of  the  static 
data  structures. 

The  sparse  LU  factorization  has  been  shown  to  be  as 
scalable  as  the  sparse  Cholesky  factorization  method. 
The  current  relea.se  of  the  SuperLU.DIST  scheme  uses 
Gaussian  elimination  with  static  pivoting,  and  the  algorithm 
is  parallelized  using  message  passing  interface  (MPl).  Other 
features  include  a  two-dimensional  (2-D)  irregular  block- 
cyclic  mapping  using  a  supemodal  structure,  and  a  kwsely 
synchronous  scheduling  with  pipelining. 

Sparse  solvers  also  deserve  mention  in  this  present  context 
because  dense  storage  of  the  coefficient  matrix  rapidly  makes 
the  solution  intractable  with  increasing  node  number.  For 
example,  a  tissue  requiring  100  x  100  x  10  nodes  would  use 
up  ~  74  GB  of  memory  for  an  8-B  floating  precision.  With  the 
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Fig.  4.  Human  head  model  from  Brooks  databa.se  with  dimensions  in 
millimeters. 


Fig.  5.  Dependence  of  runtime  for  the  144000-node  rat  model  as  a  function 
of  parallel  processors. 


present  sparse  setup,  the  actual  number  of  nonzero  entries  is 
676  000,  and  this  leads  to  less  than  6  MB. 

The  different  biomodels  currently  characterized  at  Brooks 
Base  include  rat,  frog,  monkey,  and  a  human  head.  For  compu¬ 
tational  simplicity,  we  have  chosen  to  focus  on  the  rat.  These 
models  are  based  on  anatomical  slices  from  MRI  scans.  The  rat 
model  consists  of  35  types  of  tissue,  whereas  the  human  head 
(shown  in  Fig.  4)  has  26  tissue  types.  The  spatial  resolutions  for 
the  rat  and  human  head  models  were  1  and  3  mm,  respectively. 

The  advantage  of  utilizing  parallel  processing  capability 
in  the  present  context  is  brought  out  through  the  “timing 
diagram”  in  Fig.  5.  The  figure  shows  the  factorization  time  for 
the  [M]  matrix  for  the  144C)(X)-node  rat  model  as  a  function 
of  the  processors  used  for  the  computations.  Also  shown  is  the 
average  time  taken  for  each  iterative  step  for  a  pulsed  voltage 
res|X)nse  calculation.  Two  points  are  immediately  obvious. 


Fig.  6.  Plot  of  spherical  tissue  embedded  in  a  cube  as  a  simple  example. 


First,  the  initial  factorization  process  is  very  time  consuming 
and  represents  a  large  overhead  on  the  overall  calculations. 
Any  speed-up  in  this  initial  phase  is  therefore  very  helpful. 
Employing  several  processors  greatly  reduces  this  factorization 
time.  A  monotonic  and  logarithmic  decrease  is  predicted  as 
the  number  of  processors  increa.ses  from  4  to  24.  Second,  the 
average  duration  for  each  iterative  step  of  the  time-dependent 
calculations  decreases  dramatically  as  the  number  of  processors 
is  increased  from  4  to  12.  Beyond  that,  the  average  time  for  each 
iteration  actually  increases  somewhat  and  reaches  a  saturating 
limit.  This  behavior  arises  from  a  tradeoff  between  enhanced 
average  speed  and  longer  time  spent  in  communications 
between  the  various  processors  as  their  numbers  are  increased. 
Clearly,  for  simulations  involving  short  voltage  pulses  (e.g., 
the  nanosecond  pulsars  [17]),  it  would  be  very  advantageous 
to  use  a  large  number  of  parallel  processors  as  the  factorization 
overhead  would  then  be  a  significant  percentage. 


III.  Results  and  Discussions 
A.  Test  Case  of  Uniform  Hemisphere  in  Saline 

A  test  case  was  first  analyzed  to  gauge  the  validity  of  the 
numerical  implementation  of  the  simulation  tool.  A  simple  data 
file  containing  a  single  tissue  type  (muscle)  was  chosen.  The 
geometry  was  a  spherical  mass  surrounded  by  air  in  a  cube. 
The  corresponding  3-D  plot  is  shown  in  Fig.  6.  A  slice  across 
the  bottommost  surface  (identical  to  that  at  the  other  five  faces 
of  the  cube)  extracted  from  the  data  file  is  shown  in  Fig.  7. 
The  centrally  located  circular  region  is  the  tissue.  In  the  present 
case,  electrodes  were  applied  to  the  top  and  bottom  slices  for 
the  purposes  of  obtaining  the  time-dependent  current  response. 

The  simulated  current  is  shown  in  Fig.  8  for  a  trapezoidal 
pulse  with  30-ns  rise  and  fall  times,  and  5(X)-ns  on  time. 
The  tissue  conductivity  {a)  used  was  1.33  x  10"^  {Q  •  cm)-^ 
whereas  the  dielectric  constant  was  2.7.  A  steady-state  current 
of  ~  8  A  is  predicted  during  the  pulse  on  time.  The  voltage  pro¬ 
file  across  a  cross-sectional  slice  passing  through  the  spherical 
center  and  the  two  electrodes  is  given  in  Fig.  9. 
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Fig.  7.  Plot  of  the  tissue  end-slice  for  the  data  in  Fig.  6. 


Fig.  8.  Simulated  current  response  to  trapezoidal  pulse. 


z 


Fig.  10.  2-D  projection  for  the  geometry  used  in  the  simulation. 

circular  anode  plate  is  the  solid  line  on  the  right.  The  anode 
is  parallel  to  the  cathode  and  cuts  the  sphere  in  a  circular 
disk.  In  spherical  coordinates,  the  maximum  angle  between  the 
peripheral  anode  fX)int  and  the  spherical  center  is  The 
radius  of  the  sphere  is  denoted  by  i?,  the  perpendicular  distance 
between  the  two  parallel  electrodes  is  along  the  z-axis,  and  0  is 
a  general  angle  with  <9  <  7r/2.  The  effective  resistance 
of  the  tissue  between  the  two  electrodes  can  be  computed 
analytically. 

Considering  a  differential  slice  of  thickness  dz  as  shown  in 
Fig.  10,  the  corresponding  differential  resistance  dR^fi  is 

dRefi  -  dz/[(j  Area]  (3a) 

where  the  area  of  the  differential  disk  =  7r[i?sin(^)]^.  Also. 
z  =  Rcos(0).  from  which  dz  =  -Rsui(9)d9,  and  the  total 
effective  resistance  becomes 

9  min 

RcH  =  J  j^cTTT  {i?sin(0)}^j  I .  (3b) 

-7r/2 

Canning  out  the  integration  yields  the  effective  resistance  for 
the  geometry 

=  -  {In  [tan((?min/2)]}  /{(tttR}.  (3c) 


Fig.  9.  Potential  profile  in  a  plane  through  the  electrodes  and  sphere  center. 

Next,  a  slightly  different  electrode  placement  for  the  simple 
tissue  in  saline  was  carried  out  as  a  validity  check.  The  cathode 
was  assumed  to  be  a  planar  disk  passing  through  the  center  of 
the  spherical  tissue,  with  the  anode  being  a  smaller  circular  disk 
cutting  the  sphere  parallel  to  the  cathode.  Thus,  a  2-D  projection 
of  the  geometry,  as  shown  in  Fig.  10,  resulted.  In  Fig.  10,  the 
cathode  is  represented  by  the  diameter  on  the  left,  whereas  the 


The  simulation  result  for  the  current  response  to  a  1000-V 
trapezoidal  pulse  (0.5-ns  rise  and  fall  times  and  10-ns  on  time) 
is  shown  in  Fig.  1 1 .  The  predicted  stabilized  current  during  the 
pulse  on  time  was  roughly  5 1 .5  A. 

The  effective  resistance  from  the  numerical  simulation  works 
out  to  be  1000  V/51.5  A  =  19.42  Q.  The  analytical  formula  in 
(3c)  yields  for  the  1.33  x  10“^  (D  •  cm)“^  conductivity 

fleff  =  -  {In  [tan(6lmin/2)]}  /{1.333  x  lO'®  ttIO}  n.  (4) 
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Fig.  1 1 .  Current  response  to  a  trapezoidal  applied  pulse. 


Fig.  12.  Potential  profile  at  5.775  ns  for  the  spherical  tissue. 

The  angle  0  here  is  given  by  ^  =  cos~^{7/10}  =  0.7954  rd. 
Using  this  value  in  (4)  yields  Heff  ~  20.7  This  value  is  very 
close  to  the  simulation  result  of  19.42  Q  and  roughly  validates 
the  numerical  implementation  in  this  case. 

The  voltage  profiles  for  the  above  test  simulation  are  pre¬ 
sented  in  Fig.  12  below.  The  maximum  voltage  at  the  anode  is 
1000  V  and  occurs  over  a  flat  region  intersecting  the  sphere. 
The  cathode  is  at  0  V  and  runs  through  the  parallel  diameter. 
Most  of  the  voltage  drop  occurs  between  the  two  electrodes. 
Some  electric  field  is  also  created  near  the  smaller  anode  disk. 

B.  Simulations  for  the  Rat  Model 

Next,  simulations  for  the  rat  model  are  discussed.  Fig.  13 
shows  the  bottom  view  of  a  discretized  picture  of  the  rat 
used.  The  axes  scale  are  in  millimeters,  and  the  length  of  the 
box  containing  the  rat  model  is  about  120  mm,  whereas  the 
breadth  and  depth  are  roughly  60  and  20  mm,  respectively. 
For  the  1-mm  spatial  resolution,  this  translates  into  a  total 
of  144000  nodes.  The  whole  mouse  model  contained  various 
internal  organs  and  body  parts.  Whereas  the  details  are  quite 
complicated  and  difficult  to  visualize  in  very  aspect,  a  rough 
perception  can  be  offered  by  taking  various  slices  perpendicular 
to  the  depth. 


The  contour  plots  in  Fig.  14  show  cross-sectional  slices 
at  depths  of  5,  10,  and  18  mm  from  the  top.  Corresponding 
contours  of  the  permittivity  across  the  cross  sections  are  given 
in  Fig.  15. 

The  results  of  the  internal  potential  generated  within  the  rat  in 
response  to  an  external  voltage  pulse  are  shown  and  discussed 
next.  For  simplicity,  a  trapezoidal  pulse  of  nanosecond  duration, 
with  rise  and  fall  times  of  1  ns  and  an  on  time  of  10  ns. 
was  assumed.  This  shape  and  duration  are  in  keeping  with 
the  actual  experimental  pulses  being  used  by  several  research 
groups  [31].  The  voltage  profile  at  the  midpoint  of  about 
5.775  ns  is  given  in  Fig.  16.  The  lowest  slice  with  the  four 
feet  and  tail  was  grounded.  A  uniform  potential  of  1000  V  was 
assumed  at  the  topmost  plane  parallel  to  the  ground.  As  before, 
cross-sectional  slices  were  cho.sen  perpendicular  to  the  depth 
axis  for  visualization.  A  “top  view”  showing  the  tissue  layout 
across  the  chosen  slice  is  given  in  Fig.  16,  whereas  Fig.  17 
shows  the  3-D  potential  profile  corresponding  to  the  tissue.  The 
four  “grounded”  legs  are  evident. 

A  more  realistic  pulse  with  a  waveshape  typically  generated 
by  the  pulsars  at  Brooks  City  Base  was  used  as  a  final  test  of 
the  numerical  model.  Results  for  the  applied  voltage  pulse  and 
current  are  shown  in  Fig.  18.  The  current  is  seen  to  track  the 
external  voltage  fairly  accurately.  The  displacement  currents 
are  relatively  minor.  Here,  a  comment  pertaining  to  muscular 
twitching  or  similar  electrostimulation  by  such  pulses  is  per¬ 
haps  necessary  and  relevant.  The  nerve  or  muscular  response 
is  directly  related  to  the  amount  of  charge  deposited  into  the 
biosystem  [32],  [33].  Hence,  for  maximal  effect  with  high 
efficiency,  a  monopolar  voltage  pulse  with  a  fast  rise  time  would 
appear  to  be  the  optimal  choice.  The  Brooks  pulse  shown  in 
Fig.  18  is  fairly  oscillatory  and  perhaps  not  the  most  efficient 
for  electrostimulation. 

The  time-dependent  voltage  evolution  inside  the  model  rat 
across  a  line  parallel  to  the  ground  is  shown  in  Fig.  19.  The  line 
was  located  at  a  depth  of  10  mm,  midway  along  the  longitudinal 
axis,  and  was  chosen  to  run  across  the  entire  5 1-mm  width.  The 
oscillatory  potential  at  any  point  on  the  chosen  line  is  evident. 
Although  not  shown,  the  potential  difference  across  any  chosen 
.segment  or  set  of  points  could  be  obtained  at  every  time  step. 
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Fig.  14.  Cross  sections  of  the  mouse  at  different  depths  starting  from  the  top. 
(a)  5  mm,  (b)  10  mm,  and  (c)  18  mm  along  the  depth. 

The  data  also  allowed  for  the  evaluation  of  other  parameters  of 
interest,  such  as  the  “activating  function"  for  twitching  [34],  at 
any  point  inside  the  rat. 

IV.  Conclusion 

Electrostimulation  has  been  used  for  a  variety  of  diag¬ 
nostic  and  therapeutic  applications.  Excitation  with  electrical 


Length 


Fig.  15.  Cross-scclional  contours  of  the  tissue  permittivities  for  the  rat  mcxicl. 

pulses  is  also  an  important  issue  in  the  context  of  health 
and  safety  assessment  of  various  UWB  sources.  An  important 
first  step  in  accurately  predicting  the  bioresponse  to  such  an 
electrical  stimulation  is  the  development  of  numerical  models 
for  quantifying  the  currents  and  potentials  induced  electri¬ 
cally  within  whole-body  systems.  A  side  benefit  of  such  a 
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Fig.  16.  5.775-ns  top  view  of  the  potential  profile  at  a  depth  of  12  mm. 


Fig.  17.  5.775-ns  perspective  of  the  potential  profile  at  a  depth  of  12  mm. 

simulation  is  the  capability  for  optimal  electrode  design  and 
placement. 

Typically,  whole-body  systems  are  very  complex  and  consi.st 
of  a  multitude  of  tissues,  organs,  and  subcomponents  with 
diverse  properties.  From  an  electrical  standpoint,  these  can  be 
characterized  in  terms  of  separate  conductivities  and  permittiv¬ 
ities.  Accuracy  demands  discretization  of  the  entire  simulation 
volume  into  a  fine-grained  mesh.  Practically,  this  can  place 
prohibitive  computational  requirements  on  memory  storage  and 
execution  times.  This  makes  numerical  simulations  an  arduous 
task.  Here,  we  have  presented  a  simple  yet  fast  and  efficient 
numerical  implementation  for  easing  the  computational  diffi¬ 
culty.  A  distributed  LU  decomposition  in  sparse  format  was 
employed  for  execution  on  a  cluster  of  computers  running  in 
parallel.  We  successfully  implemented  our  numerical  model 
and  successfully  tested  on  a  1-mm  spatially  resolved  rat  model. 
It  was  quantitatively  shown  that  increasing  the  processors  re¬ 
duces  the  computational  time.  Both  simple  trapezoidal  pulses 
and  more  realistic  waveshapes  from  actual  pulsars  were  used. 
Future  work  would  include  the  development  of  a  multilevel 
multigrid  capability  for  finer  detail. 


C  unen!  betv.'een  Anode  and  C  athode 


.^^plied  pulse 


Fig.  18.  Simulated  current  response  to  a  voltage  pulse  for  the  rate  biomodel 
with  the  .same  waveshape  as  typically  used  at  Brooks  City  Base. 


400  V 


Time  (ps) 

Fig.  19.  Simulated  potential  evolution  across  a  line  running  along  the  entire 
5 1-mm  width  of  the  rat  model.  The  line  was  located  at  a  depth  of  10  mm,  at  the 
midpoint  of  the  longitudinal  axis. 
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